Hessian-based sampling for high-dimensional model reduction

In this work we develop a Hessian-based sampling method for the construction of goal-oriented reduced order models with high-dimensional parameter inputs. Model reduction is known very challenging for high-dimensional parametric problems whose solutions also live in high-dimensional manifolds. However, the manifold of some quantity of interest (QoI) depending on the parametric solutions may be low-dimensional. We use the Hessian of the QoI with respect to the parameter to detect this low-dimensionality, and draw training samples by projecting the high-dimensional parameter to a low-dimensional subspace spanned by the eigenvectors of the Hessian corresponding to its dominating eigenvalues. Instead of forming the full Hessian, which is computationally intractable for a high-dimensional parameter, we employ a randomized algorithm to efficiently compute the dominating eigenpairs of the Hessian whose cost does not depend on the nominal dimension of the parameter but only on the intrinsic dimension of the QoI. We demonstrate that the Hessian-based sampling leads to much smaller errors of the reduced basis approximation for the QoI compared to a random sampling for a diffusion equation with random input obeying either uniform or Gaussian distributions

Sparse polynomial approximations for affine parametric saddle point problems

In this work we study convergence properties of sparse polynomial approximations for a class of affine parametric saddle point problems. Such problems can be found in many computational science and engineering fields, including the Stokes equations for viscous incompressible flow, mixed formulation of diffusion equations for heat conduction or groundwater flow, time-harmonic Maxwell equations for electromagnetics, etc. Due to the lack of knowledge or intrinsic randomness, the coefficients of such problems are uncertain and can often be represented or approximated by high- or countably infinite-dimensional random parameters equipped with suitable probability distributions, and the coefficients affinely depend on a series of either globally or locally supported basis functions, e.g., Karhunen--Lo\`eve expansion, piecewise polynomials, or adaptive wavelet approximations. Consequently, we are faced with solving affine parametric saddle point problems. Here we study sparse polynomial approximations of the parametric solutions, in particular sparse Taylor approximations, and their convergence properties for these parametric problems. With suitable sparsity assumptions on the parametrization, we obtain the algebraic convergence rates $O(N^{-r})$ for the sparse polynomial approximations of the parametric solutions, in cases of both globally and locally supported basis functions. We prove that $r$ depends only on a sparsity parameter in the parametrization of the random input, and in particular does not depend on the number of active parameter dimensions or the number of polynomial terms $N$. These results imply that sparse polynomial approximations can effectively break the curse of dimensionality, thereby establishing a theoretical foundation for the development and application of such practical algorithms as adaptive, least-squares, and compressive sensing constructions

On Bayesian A- and D-optimal experimental designs in infinite dimensions

We consider Bayesian linear inverse problems in infinite-dimensional separable Hilbert spaces, with a Gaussian prior measure and additive Gaussian noise model, and provide an extension of the concept of Bayesian D-optimality to the infinite-dimensional case. To this end, we derive the infinite-dimensional version of the expression for the Kullback-Leibler divergence from the posterior measure to the prior measure, which is subsequently used to derive the expression for the expected information gain. We also study the notion of Bayesian A-optimality in the infinite-dimensional setting, and extend the well known (in the finite-dimensional case) equivalence of the Bayes risk of the MAP estimator with the trace of the posterior covariance, for the Gaussian linear case, to the infinite-dimensional Hilbert space case.Comment: 16 pages, minor changes, corrected typo

Taylor approximation and variance reduction for PDE-constrained optimal control under uncertainty

In this work we develop a scalable computational framework for the solution of PDE-constrained optimal control under high-dimensional uncertainty. Specifically, we consider a mean-variance formulation of the control objective and employ a Taylor expansion with respect to the uncertain parameter either to directly approximate the control objective or as a control variate for variance reduction. The expressions for the mean and variance of the Taylor approximation are known analytically, although their evaluation requires efficient computation of the trace of the (preconditioned) Hessian of the control objective. We propose to estimate this trace by solving a generalized eigenvalue problem using a randomized algorithm that only requires the action of the Hessian on a small number of random directions. Then, the computational work does not depend on the nominal dimension of the uncertain parameter, but depends only on the effective dimension, thus ensuring scalability to high-dimensional problems. Moreover, to increase the estimation accuracy of the mean and variance of the control objective by the Taylor approximation, we use it as a control variate for variance reduction, which results in considerable computational savings (several orders of magnitude) compared to a plain Monte Carlo method. We demonstrate the accuracy, efficiency, and scalability of the proposed computational method for two examples with high-dimensional uncertain parameters: subsurface flow in a porous medium modeled as an elliptic PDE, and turbulent jet flow modeled by the Reynolds-averaged Navier--Stokes equations coupled with a nonlinear advection-diffusion equation characterizing model uncertainty. In particular, for the latter more challenging example we show scalability of our algorithm up to one million parameters resulting from discretization of the uncertain parameter field

A Nested Partitioning Scheme for Parallel Heterogeneous Clusters

Modern supercomputers are increasingly requiring the presence of accelerators and co-processors. However, it has not been easy to achieve good performance on such heterogeneous clusters. The key challenge has been to ensure good load balance and that neither the CPU nor the accelerator is left idle. Traditional approaches have offloaded entire computations to the accelerator, resulting in an idle CPU, or have opted for task-level parallelism requiring large data transfers between the CPU and the accelerator. True work-parallelism has been hard as the Accelerators cannot directly communicate with other CPUs (besides the host) and Accelerators. In this work, we present a new nested partition scheme to overcome this problem. By partitioning the work assignment on a given node asymmetrically into boundary and interior work, and assigning the interior to the accelerator, we are able to achieve excellent efficiency while ensure proper utilization of both the CPU and Accelerator resources. The problem used for evaluating the new partition is an $hp$ discontinuous Galerkin spectral element method for a coupled elastic--acoustic wave propagation problem

Weighted BFBT Preconditioner for Stokes Flow Problems with Highly Heterogeneous Viscosity

We present a weighted BFBT approximation (w-BFBT) to the inverse Schur complement of a Stokes system with highly heterogeneous viscosity. When used as part of a Schur complement-based Stokes preconditioner, we observe robust fast convergence for Stokes problems with smooth but highly varying (up to 10 orders of magnitude) viscosities, optimal algorithmic scalability with respect to mesh refinement, and only a mild dependence on the polynomial order of high-order finite element discretizations ($Q_k \times P_{k-1}^{disc}$, order $k \ge 2$). For certain difficult problems, we demonstrate numerically that w-BFBT significantly improves Stokes solver convergence over the widely used inverse viscosity-weighted pressure mass matrix approximation of the Schur complement. In addition, we derive theoretical eigenvalue bounds to prove spectral equivalence of w-BFBT. Using detailed numerical experiments, we discuss modifications to w-BFBT at Dirichlet boundaries that decrease the number of iterations. The overall algorithmic performance of the Stokes solver is governed by the efficacy of w-BFBT as a Schur complement approximation and, in addition, by our parallel hybrid spectral-geometric-algebraic multigrid (HMG) method, which we use to approximate the inverses of the viscous block and variable-coefficient pressure Poisson operators within w-BFBT. Building on the scalability of HMG, our Stokes solver achieves a parallel efficiency of 90% while weak scaling over a more than 600-fold increase from 48 to all 30,000 cores of TACC's Lonestar 5 supercomputer.Comment: To appear in SIAM Journal on Scientific Computin

Hessian-based adaptive sparse quadrature for infinite-dimensional Bayesian inverse problems

In this work we propose and analyze a Hessian-based adaptive sparse quadrature to compute infinite-dimensional integrals with respect to the posterior distribution in the context of Bayesian inverse problems with Gaussian prior. Due to the concentration of the posterior distribution in the domain of the prior distribution, a prior-based parametrization and sparse quadrature may fail to capture the posterior distribution and lead to erroneous evaluation results. By using a parametrization based on the Hessian of the negative log-posterior, the adaptive sparse quadrature can effectively allocate the quadrature points according to the posterior distribution. A dimension-independent convergence rate of the proposed method is established under certain assumptions on the Gaussian prior and the integrands. Dimension-independent and faster convergence than $O(N^{-1/2})$ is demonstrated for a linear as well as a nonlinear inverse problem whose posterior distribution can be effectively approximated by a Gaussian distribution at the MAP point

Solution of nonlinear Stokes equations discretized by high-order finite elements on nonconforming and anisotropic meshes, with application to ice sheet dynamics

Motivated by the need for efficient and accurate simulation of the dynamics of the polar ice sheets, we design high-order finite element discretizations and scalable solvers for the solution of nonlinear incompressible Stokes equations. We focus on power-law, shear thinning rheologies used in modeling ice dynamics and other geophysical flows. We use nonconforming hexahedral meshes and the conforming inf-sup stable finite element velocity-pressure pairings $\mathbb{Q}_k\times \mathbb{Q}^\text{disc}_{k-2}$ or $\mathbb{Q}_k \times \mathbb{P}^\text{disc}_{k-1}$. To solve the nonlinear equations, we propose a Newton-Krylov method with a block upper triangular preconditioner for the linearized Stokes systems. The diagonal blocks of this preconditioner are sparse approximations of the (1,1)-block and of its Schur complement. The (1,1)-block is approximated using linear finite elements based on the nodes of the high-order discretization, and the application of its inverse is approximated using algebraic multigrid with an incomplete factorization smoother. This preconditioner is designed to be efficient on anisotropic meshes, which are necessary to match the high aspect ratio domains typical for ice sheets. We develop and make available extensions to two libraries---a hybrid meshing scheme for the p4est parallel AMR library, and a modified smoothed aggregation scheme for PETSc---to improve their support for solving PDEs in high aspect ratio domains. In a numerical study, we find that our solver yields fast convergence that is independent of the element aspect ratio, the occurrence of nonconforming interfaces, and of mesh refinement, and that depends only weakly on the polynomial finite element order. We simulate the ice flow in a realistic description of the Antarctic ice sheet derived from field data, and study the parallel scalability of our solver for problems with up to 383M unknowns.Comment: 31 page

A randomized maximum a posterior method for posterior sampling of high dimensional nonlinear Bayesian inverse problems

We present a randomized maximum a posteriori (rMAP) method for generating approximate samples of posteriors in high dimensional Bayesian inverse problems governed by large-scale forward problems. We derive the rMAP approach by: 1) casting the problem of computing the MAP point as a stochastic optimization problem; 2) interchanging optimization and expectation; and 3) approximating the expectation with a Monte Carlo method. For a specific randomized data and prior mean, rMAP reduces to the maximum likelihood approach (RML). It can also be viewed as an iterative stochastic Newton method. An analysis of the convergence of the rMAP samples is carried out for both linear and nonlinear inverse problems. Each rMAP sample requires solution of a PDE-constrained optimization problem; to solve these problems, we employ a state-of-the-art trust region inexact Newton conjugate gradient method with sensitivity-based warm starts. An approximate Metropolization approach is presented to reduce the bias in rMAP samples. Various numerical methods will be presented to demonstrate the potential of the rMAP approach in posterior sampling of nonlinear Bayesian inverse problems in high dimensions

Inexact Newton Methods for Stochastic Nonconvex Optimization with Applications to Neural Network Training

We study stochastic inexact Newton methods and consider their application in nonconvex settings. Building on the work of [R. Bollapragada, R. H. Byrd, and J. Nocedal, IMA Journal of Numerical Analysis, 39 (2018), pp. 545--578] we derive bounds for convergence rates in expected value for stochastic low rank Newton methods, and stochastic inexact Newton Krylov methods. These bounds quantify the errors incurred in subsampling the Hessian and gradient, as well as in approximating the Newton linear solve, and in choosing regularization and step length parameters. We deploy these methods in training convolutional autoencoders for the MNIST and CIFAR10 data sets. Numerical results demonstrate that, relative to first order methods, these stochastic inexact Newton methods often converge faster, are more cost-effective, and generalize better